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Abstract 



High-resolution simulations within the GOY shell model are used to study various 
scaling relations for turbulence. A power-law relation between the second-order in- 
termittency correction and the crossover from the inertial to the dissipation range 
is confirmed. Evidence is found for the intermediate viscous dissipation range pro- 
posed by Frisch and Vergassola. It is emphasized that insufficient dissipation-range 
resolution systematically drives the energy spectrum towards statistical-mechanical 
equipartition. In fully resolved simulations the inertial-range scaling exponents de- 
pend on both model parameters; in particular, there is no evidence that the conser- 
vation of a helicity-like quantity leads to universal exponents. 
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1 Introduction 

One of the most challenging aspects of turbulent fluid flow is the rapid rise 
in the number of active degrees of freedom as the Reynolds number increases. 
Given present day computational resources, this renders direct numerical sim- 
ulation of the asymptotic high-Reynolds-number regime of fully developed 
turbulence intractable. Various reduced models of the dynamics have there- 
fore been developed to capture the essential features of the turbulent cascade 
of energy from large to small scales [22,47-49,40,46,16,17,10,4]. These models 
reflect many of the statistical features expected to influence the energy trans- 
fer in realistic flows [46,23-25]. They have been used to study intermittency, 
multifractal behavior, and the scaling of energy and dissipation (see also [7,8]). 
They have also been extended to include passive scalars [30], magnetic fields 
[6], and even polymers [2]. 

In this note we focus on the GOY shell model, so named after its initial pro- 
ponents Gledzer [22] and Ohkitani and Yamada [47-49] (see also [10,4] for 
reviews). We consider within high- resolution, fully resolved, and converged 
numerical simulations four related issues: the Reynolds number dependence of 
the crossover to the dissipation scale; the existence of an intermediate dissi- 
pation range; the effect of modal truncation in the viscous range on inertial- 
range scaling; and the relation between inertial-range scaling and the presence 
of a second conserved helicity-like quantity. Some of these aspects have been 
discussed for shell models in isolation before (e.g. [35,19,36,37]). Key observa- 
tions of the present paper are: the crossover from inertial to dissipative scales 
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takes place at a length much smaller than the Kolmogorov estimate; these 
scales must be properly resolved to avoid contaminating the inertial range; 
and, finally, the high-resolution, fully converged simulations presented here do 
not support the conjecture [31] that inertial-range exponents in the helicity- 
preserving GOY model are universal. 

The usual estimate for the scale on which dissipation begins to dominate, 
based on the energy dissipation density e and the kinematic viscosity u, is the 
Kolmogorov scale 



However, numerical simulations of passive scalars show that in order to recover 
small-scale fluctuations, gradients, and dissipation accurately, one has to go 
well below this scale [43]. A straightforward line of reasoning expanded on in 
Sect. 3 shows that in the presence of intermittency corrections, as reflected 
in the energy density scaling like k~ 5 ^ 3 ~ s with < 5 < 4/3, the dissipative 
range sets in at a scale rj d that is a fraction of rj that decreases as the Reynolds 
number increases (i.e. as the viscosity decreases): 



The integral or outer scale L appearing here is the characteristic scale of 
the external stirring force. A more detailed argument for the variation of 
the energy density in the transition between the inertial and viscous ranges, 
based on the multifractal model of intermittency, shows that an intermediate 
dissipation range develops that widens with increasing Reynolds number [21]. 
In Sect. 4, we present high-resolution numerical evidence for the existence of 
an intermediate dissipation range in the GOY model. However, this part of 
the spectrum is so steep that it contributes negligibly to the energy balances 
used in Sect. 3 to determine the dissipation wavenumber. 

As part of these studies we also considered the influence of modal truncations 
on spectra. If the modes are truncated very far in the viscous range, it is 
reasonable to expect that the behavior will not change. However, as the cutoff 
moves closer and closer to the dissipation scale, it begins to inhibit sufficient 
dissipation. As a consequence, the pileup of energy close to the transition 
between the inertial- and viscous-dominated regimes (the bottleneck effect, 
which also occurs in simulations of the Navier-Stokes equation [19,36,37]) 
is enhanced to the point where the system tends towards a weakly damped 
equilibrium: the distribution of energy approaches a 1/k scaling, corresponding 
to an equipartition of the energy contents of the geometrically spaced shells. 
A similar kind of effect has been noted in models with hyperviscosity, where 
the rapid quenching of high-wave number modes also effects the inertial-range 
scaling [35]. 
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In Sect. 5, using fully resolved simulations, we revisit the studies of Kadanoff 
et al. [31,42]. These authors considered different parameters for the GOY 
model and reported that the intermittency corrections fall into two categories, 
depending on whether or not the model preserves a helicity-like quadratic 
quantity. On revisiting these studies with higher-resolution simulations, we 
do not find such a correspondence, but rather a persistent dependence of the 
intermittency corrections on all parameters of the model. 

Before addressing these issues, it is helpful to describe in the next section the 
underlying model and some related technical details. 



2 Preliminaries 



The GOY model has, as primary dynamical variables, complex velocities u n 
representing the velocity field in shell number n, where n = 0, 1, . . . , N — 1. 
They evolve according to 




<S n + F n , (3) 



where 



S n = ik n (au* n+1 u*n+2 + jK-l u *n+l + J<-l<- 2 j , ( 4 ) 

given the boundary conditions w_2 = «-i = un = %+i = 0. The wavenum- 
bers k n = /c A n scale geometrically. It is customary to rescale time so that 
a — 1 and require that a+P+j = 0, so that the nonlinear terms in S n conserve 
the energy | YJ n \u n \ 2 . A second invariant | YJ n A;^|-u„| 2 is also conserved, where 
p = — log A (— (3 — 1). Of particular interest is the case where A = 1/(1 + (3), 
when this invariant takes the form of a quantity H — | YJ n (— l) n k n \u n \ 2 with 
the same dimensions and sign indefiniteness as the helicity invariant of three- 
dimensional Navier Stokes turbulence. 

For the standard case considered by Kadanoff et al. [31,42], ko = 2~ 3 , A = 2, 
f3 — 7 = — l/2, v = 10~ 7 , N = 22 and the forcing is confined to shell 3: 

F n = fS n , 3 . (5) 

The amplitude of the external forcing is constant: / = 5(1 + i) x 10~ 3 . In 
these models, one observes an inertial-range power-law scaling of the mean 
shell energy \(\u n \ 2 ) ~ k~ 2 ^ 3 , corresponding to the energy spectrum E(k n ) = 
\{\ u n\ 2 ) I (k n +i~k n ) ~ ^« 5//3 and reminiscent of the "K41" scaling predicted by 
the Kolmogorov theory [32]. Here () denotes an ensemble average. At higher 
shell indices one observes an exponential decrease of the shell amplitudes, 
corresponding to a viscous range [1]. 
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On multiplying (3) by u* n and summing over shells m to N — 1, one arrives at 
the energy balance 

d N ~ 1 

-J, \( u n)\ 2 = n m - e m , (6) 

n=m 

where IT m = 2 Re J2n=m(^nUn) represents the transfer of energy by the non- 
linear terms into shells m to N — 1 and 

e m = 2v X>n<K| 2 } -2 Re £<FX>, (7) 

is the total transfer, m« dissipation and forcing, out of this range of shells. A 
positive value for U m represents a flow of energy to shells m and higher. When 
v — and F = 0, energy conservation implies that 

m— 1 

n» = -2E(W. ( 8 ) 

n=0 

so that n = = 0. In a statistical steady state, the left-hand side of (6) 
vanishes and Il m = e m . This condition serves as an excellent numerical di- 
agnostic for discerning when a steady state has been reached. In practice, 
appealing to the ergodic theorem, running time averages (say from T/2 to the 
final time T) are used to approximate the required ensemble averages. 

As pointed out in [41], it is highly advantageous to use an exponential in- 
tegrator to solve (3) exactly on the linear time scale. Exponential integra- 
tors [13,28,3,15,29] are well suited to linearly stiff systems that exhibit dy- 
namics on a wide range of time scales. They allow for much larger time steps 
than either classical nonstiff methods or simple integrator factor (Rosenbrock) 
methods [26] . However, instead of using the second-order exponential Adams- 
Bashforth scheme used by Pisarenko et al. [41] or the fixed time step schemes 
of Cox and Matthews [15], we used the exponential version [11] of the adap- 
tive third-order Bogacki-Shampine [9] Runge-Kutta integrator described in 
Appendix A. With a nonstochastic forcing, this scheme can reuse the final 
source evaluation from the previous time step. 

Instead of a constant forcing, it is also possible to apply a central white- 
noise stochastic forcing, where / is delta-correlated in time. The advantage of 
this forcing is that it yields a prescribed mean energy injection e = |(|/| 2 ), 
independent of the lattice size and forcing [39]. 

The fully resolved spectra and cumulative energy transfers Il n and e n corre- 
sponding to a series of decreasing values of v are shown in Fig. 1. Here we 
set ko — 1 and used the standard values (/3, A) = (—1/2,2) with a complex- 
valued white-noise random forcing restricted to the first shell, such that the 
total energy injection e is 1. Initially all shells were assigned the same nonzero 
real amplitude (a statistical mechanical equipartition of energy); the complex 
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forcing and nonlinearity then cause them to develop imaginary parts. The 
dashed vertical lines in Fig. 1 separate regions of equal energy dissipation, 
providing a convenient definition for the dissipation wavenumber k^, namely, 
2v$ d k 2 E(k) dk = e/2. As seen by the near coincidence of with the maxima 
of the energy dissipation, this definition of kd closely approximates a point of 
inflection, with respect to shell number, of the cumulative energy transfer [34]. 

The highest resolution case was run 1.2 x 10 9 adaptive time steps for a total 
of 9756 time units, or roughly 4000 large-eddy turnover times. In Fig. 2, we 
illustrate the power-law behavior of kd vs. v evident in Fig. 1; using a least- 
squares fit, we find that 

kd ~ i/" " 7856 , (9) 

with the exponent determined to within a statistical error of 0.0005. Thus, 
instead of the K41 value —0.75 arising from (1), we obtain an exponent much 
closer to the anomalous value —0.775 deduced from (2), using the second-order 
intermittency correction S = 0.0438 measured in section 3. 

The slow decrease with wavenumber of the energy spectra compensated by the 
simple algebraic scaling indicates an intermittency correction to the exponent. 
The relationship between this intermittency correction and the dissipation 
wavenumber kd will be considered in the next section. 



3 Intermittency-length scale paradox 



Deviations of the p-th order structure function scaling exponents from the 
Kolmogorov value p/3 require the presence of a length scale £ for dimensional 
correctness: 

(kjy 6p ~ e*/V^~*» ( 10 ) 

where ( p = p/3 + 5 p . Two natural choices for this length scale are the external, 
stirring scale L and the small scale r] d at which dissipation terminates the 
inertial range. Concerning this choice, we now give two elementary arguments 
that at first glance lead to seemingly conflicting results. 

The energy spectrum may be expressed in terms of the intermittency exponent 
5 = 5 2 : 

E{k) = Ce 2/3 k- 5 / 3 (k£)- 5 . (11) 

On balancing the energy injection e and energy dissipation between the largest 
scale L and some dissipation scale rjd, one finds 
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Fig. 1. (a) Compensated energy densities for the GOY model driven by a white-noise 
random forcing, normalized by the K41 law, for v = 10 -3 , 10 -4 , 10 -5 , 10 -6 , 10 -7 , 
10" 8 , 1(T 9 , 10" 10 , and 10 -11 . (b) Cumulative energy transfer functions: the dashed 
(dotted) curves indicate the value of Il n (e n ) at the wavenumbers k n . The energy 
dissipation z/fc 2 (|« n | 2 ) (solid curves) is maximal very near kd (vertical dashed lines). 



= 2u k 2 E{k) dk = 2Cue 2/3 r 5 f 2 *'™ k 1/3 - s dk 

J2n/L J2n/L 
( 1 \ 4/3 " 5 

~ ' IW (^« L )' ( 12 ) 



where we have tacitly assumed that S < 4/3, so that the upper integration 
limit dominates. If one takes r\d to be the Kolmogorov scale r\ = (z/ 3 /e) 1//4 , 
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Fig. 2. Scaling of kd vs. v for the simulations in Fig. 1. The line joining the data 
points is a least-squares fit. 

then on comparing the left- and right-hand sides of (12), one obtains 



For 5^0, such a balance is possible only if i ~ 77; that is, the relevant choice 
for t is the dissipation scale. 

On the other hand, the Cauchy-Schwarz inequality implies that structure 
functions satisfy [20] 



So if I is the smallest excited length scale {£ ~ rj), then r/£ > 1 for all r, which 
would imply that S p (and hence ( p ) is convex. Alternatively, if £ is the largest 
excited length scale (£ ~ L), then r/£ < 1 for all r, and 5 P (and hence ( p ) is 
concave. 

Both Navier-Stokes turbulence and the GOY model have concave exponents, 
as illustrated in Fig. 3 and tabulated in Table 1 for (/3 = —1/2, A = 2). 
In fact, the curvature of ( p is very close to that observed experimentally for 
three-dimensional Navier-Stokes turbulence [27]. This supports i = L (see 
also Ref. [18]), even though (13) appears to require that £ = r\. The above 




(13) 




or 




(16) 



1 2 3 4 5 6 

P 

Fig. 3. Concavity of structure functions exponents C, p vs. p for (3 = —0.5 and 
A = 2. The circles represent numerical values obtained by a least-squares fit of 
the flux (21) between k = 10 2 and 10 6 for the case v = 10 . The dashed line 
is the Kolmogorov prediction £ p = p/3; the solid line indicates the She-Leveque 
scaling ( p = p/9 + 2 - 2(2/3)p/ 3 , and the squares represent the experimental values 
measured for 3D turbulence by van de Water & Herweijer [27]. In the inset we graph 
the dimension D(h) determined from our numerical £ p values by inversion of the 
multifractal Legendre transformation. 
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2 


3 


4 


5 


6 




0.7105 


0.9999 


1.262 


1.503 


1.726 




±0.0005 


±0.0005 


±0.001 


±0.001 


±0.002 



Table 1 

Numerically computed structure function exponents £ p and the statistical errors A^ p 
(from a least-squares fit between k = 10 2 and 10 6 ) for f3 = -0.5 and A = 2. 

calculations assume spectra with a sharp cut-off at the dissipation scale, but 
the same analysis with the same conclusion can be performed for spectra with 
an exponential cut-off. 

This apparent contradiction is easily resolved: a steeper-than-Kolmogorov 
spectrum necessitates integrating to scales smaller than 77 to obtain sufficient 
dissipation [21]. If we set r] d ~ 77 {^j , we obtain the energy balance 



/ Tj \ 4/3) 



(17) 



This is consistent with 



L if S + a(S - 4/3) = 0, that is, if 
1. 



a 
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For < 5 < 4/3, we see that a > 0. The conclusion is that the inertial range 
extends to a scale i] d that is smaller than the Kolmogorov scale r\: 



(jj-il" (y) 4 , (19) 



for a fixed outer scale L. When the energy injection e is also fixed, we thus 
have 



Vd 



u 1 '^. (20) 



Equivalently, on setting t = L, one can obtain that this scaling directly 
from (12). 

To verify this relation numerically, we followed Kadanoff et al. [42] and com- 
puted ( p via a least-squares fit of the time-averaged pth-order flux 



y 



T ( 

lm I u n u n+ iu n+2 H — u n -iu n u n+ i 



P/3\ 

>, (21) 



in order to filter out the well-known period-three oscillation in solutions of the 
GOY model [41]. We chose the interval [10 2 , 10 6 ] for doing the least-squares 
fit, as the logarithmic slopes of E n>p were nearly flat in this region and, in the 
case n — 3, very close to the exact value of —1. From this fit, we obtained 
the value C2 = 0.7105 ± 0.0005 tabulated in Table 1. The error quoted here 
and elsewhere is the statistical error from the least-squares fit to the data. 
Other significant sources of error are the choice of wavenumber interval for 
the fit and the time- averaging window. By varying the fitting interval and 
time-averaging window, we confirmed that the error from these sources has 
roughly the same magnitude as the statistical fitting errors (cf. Fig. 12). On 
substituting the measured second-order intermittency correction 82 = (2 — 
2/3 = 0.0438 into (20), we obtain 

k d - u- - 775 , (22) 

in remarkable agreement with the observed dissipation wavenumber scaling (9). 
This result has implications for theoretical predictions of how the dissipation 
wavenumber should depend on Reynolds number. Practically, it suggests for 
numerical simulations how the maximum truncation wavenumber should scale 
as the viscosity is decreased. Our results in section 5 indicate, for the standard 
case of Ref. [31], that the upper truncation wavenumber must be chosen to 
be at least three times larger than the dissipation wavenumber k d in order to 
obtain reliable inertial-range dynamics. 
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4 Intermediate dissipation range 



Fluctuations in the energy transport and dissipation were previously consid- 
ered by Frisch and Vergassola within a multifractal model [21] and studied 
numerically by Nakayama [38]. Frisch and Vergassola noted that the inertial- 
range scaling exponent h for velocity differences has an influence on the 
spectrum at the dissipative scale as well. Equating the viscous time and the 
turnover time, they found an /i-dependent viscous cut-off, 

V (h) ~ ^/(i+ft). (23) 

This agrees with our equation (20) after the identification h — 1/3 — 5. Frisch 
and Vergassola appealed to a geometric characterization of the sets of points 
that contribute to a certain velocity scaling exponent h for a given separation 
scale and introduced the Hausdorff dimension D(h) of these sets. In the present 
shell model, such an interpretation is not available; nevertheless, the spectrum 
can still be expressed in terms of an effective dimension 

D(h) = inf (p/i + 3 - Q (24) 

calculated directly from the scaling exponents ( p , as displayed in the inset of 
Fig. 3. The dominating scaling exponent h is given by the slope of the graph 
of C p at p = 2 in Fig. 3. On fitting a Bezier cubic spline through the numerical 
data points, we find that h = 0.309, which yields 

kd ~ v' ™. (25) 

In comparison with the measured scaling (9), this result is slightly less accurate 
than (22). If we use the She-Leveque expression [45] for ( p instead of our 
measured values, we find h = 0.317, which yields the dissipation wavenumber 
scaling exponent —0.759. 

We show in Fig. 4 that the rescaling function proposed by Frisch and Vergas- 
sola collapses all of the spectra in Fig. 1 to a single curve. In Fig. 5, we use 
our computed function D(h) to provide explicit numerical evidence for the 
existence of intermediate dissipation ranges in these simulations. We compare 
the portion of each spectrum in Fig. 1 for k > 1/rj with the intermediate 
dissipation range spectrum proposed by Frisch and Vergassola [21]: 

E(k) = 4jfe-4-2fc(k)+fl(h(k)) ( fc > 1 / r? ) ) (26) 

where h(k) = — 1— log vj log A;. We set the constant of proportionality A to 0.5. 
For the lower resolution runs, one notices only a narrow intermediate range, 
with the numerical spectra quickly entering the full dissipation range. At the 
highest resolutions, the intermediate dissipation range appears to consist of all 
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resolved wavenumbers above 1/rj. Since the energy in these intermediate dis- 
sipation ranges decays very rapidly beyond ka, it contributes negligibly to the 
energy balance (17). In Fig. 6, we illustrate the rapid growth of intermittency 
in the intermediate dissipation range by plotting the flux kurtosis E n 4 /E^ 2 
vs. A n , consistent with the findings of Ref. [14], although a direct comparison 
of our deterministic model with the random cascade model studied there is 
not available. Similar (but less-smooth, due to the period-three oscillation) 
results were obtained for the kurtosis of the shell velocities u„. 



bo 

o 



o 




0.4 0.6 0, 

logfc/log(l/i/) 



Fig. 4. Collapse of the spectra in Fig. 1 with the rescaling proposed by Frisch and 
Vergassola. 



CO 
ITS 




Fig. 5. Intermediate dissipation ranges in the simulations of Fig. 1. For each case, 
the dots indicate the portion of the numerical spectrum for k > 1/rj, the solid line 
depicts the spectrum given by (26), and the vertical dashed line separates regions 
of equal energy dissipation. 
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Fig. 6. The rapid growth of the flux kurtosis for the simulations of Fig. 1 in each 
intermediate dissipation range, beginning at kj (denoted by a vertical dashed line). 

In summary, we were able to determine the effective dimension D(h) and 
study its effects on the dissipation-range energy spectrum. The steepening of 
the energy spectrum in the inertial range implies that the dissipation range 
sets in well after the Kolmogorov scale is reached, namely at the scale deter- 
mined by (20) or (23). Frisch and Vergassola [21] note that this reduction in 
scale compared to the Kolmogorov scale is relatively small, but our numerical 
simulations are precise enough to resolve it. 



5 Dissipation-scale effects on inertial-range dynamics 

When the number of shells is kept fixed, and the viscosity is reduced, the 
viscous range will become narrower and narrower, until it finally disappears. 
Schorghofer et al. [42] call this the regime of strong interaction. However, as 
this happens, the effect of the lack of resolution no longer remains confined to 
the small scales: it begins to affect the inertial-range intermittency corrections 
as well. This effect is illustrated in Fig. 7 for simulations with fixed viscosity 
and varying N. In particular, we see on comparing the N = 18 and N = 20 
spectra that it is necessary to resolve the small-scale dynamics well beyond the 
dissipation wavenumber (denoted by the vertical dashed line). The smallest 
number of shells for which the large-scale dynamics was adequately resolved 
was N = 20, which requires a maximum wavenumber more than three times 
higher than k&. For the case iV = 14, where the viscous range was completely 
resolved, it was necessary to use the conservative C-RK5 integrator described 
in Appendix B to discretize the nonlinear terms in a manner that respects 
exact energy conservation. Otherwise, as illustrated for the nonconservative 
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solution in the transfer function inset, viscous dissipation will be insufficient 
to absorb the temporal energy discretization error, so that a global energy 
balance can never be reached. In Fig. 8, we use this conservative integrator 
to illustrate the eventual equipartition of shells energies and resulting k~ x 
energy spectrum, in the extreme limit of zero dissipation, given an initial k~ 2 
spectrum. 




10" 1 10° 10 1 10 2 10 3 10 4 

k 



10 5 10 6 



Fig. 7. The effect of a finite number of shells on the spectrum for the standard case 
(P = -1/2, A = 2) of Ref. [31]. For N = 26 (dashed) and N = 20 (dotted), there 
is no discernible difference between the energy spectra, showing that a sufficient 
number of shells in the viscous subrange is included. The vertical dashed line denotes 
the dissipation wavenumber kj- For N = 18 (solid), the dissipation range is barely 
resolved and the inertial range is contaminated. For N = 14, the inertial range tends 
toward the statistical equipartition spectrum k~ l : the short solid curve indicates the 
spectrum obtained with a conservative integrator; the long-dashed curve indicates 
the spectrum obtained with a nonconservative exponential integrator that violates 
global energy balance. 

In view of these simulations, the conjecture in [42] on the role of preserva- 
tion of the helicity H has to be reconsidered. Kadanoff et al. noted that their 
simulations for the standard case (/3 = — 1/2, A = 2) closely reproduced exper- 
imentally measured intermittency exponents of three-dimensional turbulence 
(cf. Fig. 3). At the time, it had already been observed [5] that for fixed A, the 
intermittency corrections strongly depend on the parameter (3, as illustrated 
in Fig. 9. Kadanoff et al. then hypothesized that experimentally observed ex- 
ponents should be obtained on the curve A = 1/(1-1-/3) along which the second 
conserved invariant of the GOY model is analogous to the helicity invariant 
of three-dimensional Navier-Stokes turbulence. 



However, on repeating their simulation with a sufficient number of shells, 
(particularly the case with (3 = —3/10 and A = 10/7, which requires 44 rather 
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Fig. 8. Equipartition spectrum for 16 shells obtained with the conservative C-RK5 
integrator with e = v = 0. 




Fig. 9. Scaling anomalies S p , with A = 2, for p = 2, 3, 4, 5, and 6 (from top to 
bottom) vs. f3, as obtained by a least-squares fit between k = 500 and 2 x 10 6 . The 
number of shells is 32 and the viscosity v is 10 -11 . The system is driven with a 
stochastic white-noise force on the first shell such that the energy injection e = 1. 

than 22 shells), we find that the structure function exponents change dramat- 
ically. In particular, in view of Fig. 10, it appears that the distinction made 
between the cases with and without helicity conservation in Fig. 4 of Ref. [31] 
cannot be upheld: helicity preservation does not lead to unique intermittency 
corrections. This is further emphasized in Fig. 11, where the variation of the 
anomalies S p with f3 along the helicity preserving curve is shown. Fig. 12 em- 
phasizes that the computed slopes are well resolved, and independent of the 
fitting interval, particularly for small values of j3 (corresponding to small A 
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0.04 


0.05 
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0.043 


0.044 


0.042 




-0.18 


-0.11 


-0.08 


-0.07 


-0.08 


-0.07 




-0.75 


-0.50 


-0.30 


-0.27 


-0.32 


-0.30 



Table 2 

Values of /?, A, N, 62, 64, and 5q for the helicity-preserving cases. 

and large N), where the departure from the conjecture is seen to be most 
significant. 

While Fig. 10 might appear to support the conjecture at least for p = 2, we 
see that neither Fig. 11 nor the higher-resolution analog Fig. 13 of Fig. 10 
generated using data from Figs. 9 and 11, agree with it. Our high- resolution 
numerical data set is sufficiently well resolved that we did not need to invoke 
extended self-similarity to obtain the results in Fig. 13. 

0.5 



0.4 



0.3 



0.1 
0.0 
-0.1" 

2 3 4 5 6 

P 

Fig. 10. Recalculation of Fig. 4 from Kadanoff et al, fitting between k = 20 and 
200. The first three (— /3, A) pairs respect helicity conservation; the fourth does not. 
All four cases exhibit large anomalies, independent of whether helicity is preserved. 



6 Conclusions 

In this work, we exploited the numerical advantages of the GOY model, com- 
bined with efficient integration algorithms, to obtain fully converged spectra 
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Fig. 11. Scaling anomalies 5 P from top to bottom for p = 2, 3, 4, 5, and 6 in the 
helicity-preserving case where A = 1/(1 + (3) for the values of (3 and number of shells 
N listed in Table 2, as obtained by a least-squares fit between k = 500 and 2 x 10 6 . 
The viscosity v is 10 -11 . The system is driven with a stochastic white-noise force 
on the first shell such that the energy injection e = 1. 
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Fig. 12. Determination of 5q from the logarithmic slope of E nj 6 for different values 
of p, with A = 1/(1 + 13). 

that allowed us to address links between dissipation, scaling, and conserved 
quantities. 

We documented that the typical dissipation scale is much smaller than the 
Kolmogorov length, due to the inherent steepening of the spectrum by inter- 
mittency corrections. While the detailed form of the spectrum beyond this 
point is influenced by the probability distribution of the velocity scaling ex- 
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2 3 4 5 6 

P 

Fig. 13. Analog of Fig. 10 generated from the high-resolution data in Figs. 9 and 11. 

ponent, as discussed by Frisch and Vergassola [21], the delayed transition is 
not. Overall, the effect is relatively small, but easily identified in simulations 
of the GOY model. 

The shell model simulations presented here also underscore the need to fully 
resolve the small scales. They demonstrate how insufficient resolution of the 
viscous subrange can affect inertial-range scaling exponents. In particular, 
when the viscous range is fully resolved, one observes significant variations in 
the spectral exponents even along the helicity preserving curve, thus disproving 
a previous conjecture [31]. For the standard case of Ref. [31], we found that 
the upper truncation wavenumber must be at least three times higher than 
the dissipation wavenumber kd in order to capture the proper inertial-range 
dynamics. 

While current direct numerical simulations of three-dimensional Navier-Stokes 
turbulence barely reach into the inertial range, the implication of this rela- 
tion between the second-order intermittency correction and the dissipation 
wavenumber is that, for future numerical simulations, it will become increas- 
ingly important to account properly for unresolved viscous scales. 
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A Exponential adaptive (3,2) Bogacki-Shampine Runge-Kutta in- 
tegrator 



The GOY model can be written as a set of ordinary differential equations for 
the components of the vector u = (u , u±, . . . , Wjv-i) : 

— + vu = S(u), (A.l) 

where v is a diagonal matrix with entries (iz/cq, vkl, . . . , vk%_^). 

A general s-stage autonomous Runge-Kutta scheme that evolves an initial 
value uq by a time step r to a final value u s may be written in the form 

j-i 

Ui = u + rJ2 a ij S ( u j) (i = l,...,s). (A. 2) 

3=0 



Letting x = —tv, the 4-stage exponential Bogacki-Shampine Runge-Kutta 
order (3,2) embedded pair [11] uses the coefficients 



010 =2^ \2 X 



i 



9 /3 \ 3 /l 



a 20 = -g>\ j - a 21 , 021 = -v? 2 j + -if 2 y-x 

1 4 2 

a 30 =<pi(x) - a 31 - a 32 , a 3 i = ^(pi{x), a 32 = -(p 2 {x) - -<Pi(x), 

17 1 2 1 

a4o = (pi(x) - —<f 2 (x), o 4 i = -<f 2 (x), a i2 = -ip 2 (x), a 43 = -<f 2 (x), 

where the (here diagonal) matrix functions tpi and (p 2 are given by tpi(x) = 
x~ l (e x — 1) and tp 2 (x) = x~ 2 (e x — 1 — x). The value n 3 is the third-order 
solution, while the value u 4 provides a second-order estimate that can be used 
to optimally adjust the time step. Since S(u 3 ) is just the value of S at the 
initial u for the next time step, no additional source evaluation is actually 
required to compute (unless an extra stochastic forcing term is added to 
u 3 before proceeding to the next time step). 

Care must be exercised when evaluating the diagonal components (pi(x) and 
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(fi 2 (x) near 0; optimized double precision routines for evaluating these functions 
are available at www. math. ualberta. ca/~bowman/phi .h. 



B Conservative adaptive (5,4) Runge-Kutta integrator 



Here we describe a fifth-order conservative integrator (C-RK5) similar to the 
second-order conservative predictor-corrector algorithm introduced in [12,44,33]. 
We first consider the case where each component of u in (A.l) is real. The 
scheme shares the same first five stages as the classical 6-stage Cash-Karp 
Fehlberg Runge-Kutta (4,5) embedded (RK5) pair. However, the final conser- 
vative stage is derived by applying the final conventional stage to the evolution 
equation for each (say, r th ) component u of u: 



d^_ 

dt 



2S(u)u, 



(B.l) 



where S(u) = S r (u)—v r u r . One then transforms the resulting estimate for the 
new value of u 2 back into an estimate for u. The conventional RK5 estimate 
can be used to resolve the appropriate branch of the square root: 



u 6 



Sgn (u + T a 5jS(Uj) 



3=0 



\ 



2t jo 5j %; 



u 



3 ' 



(B.2) 



3=0 



In regions where the argument of the square root becomes negative, the time 
step needs to be temporarily reduced. The resulting solution is then compared 
with the conventional fourth-order estimate uq + TY^=oO,' 5 jS(uj) (for a sec- 
ond set of coefficients a' 5 j) to adjust the time step appropriately. When u is 
complex, one applies (B.2) separately to the real and imaginary parts of u. 

The transfer function Ii M is based on partial sums of the time averages of S n u* n 
from time t 1 to t 2 . We compute these time averages as [T n (t 2 )—T n (tx)]/ (t 2 —ti), 
where 



dT n 
dt 



= S n u* n . 



(B.3) 



Energy conservation relies on the fact that J2n=o T n is invariant (since Y,n=o SnU* n 
0). Since integrators of the form (A. 2) conserve such invariants that are linear 
in the evolved variable [33], one can use a conventional RK5 integrator to 
conservatively evolve (B.3). 



We compute time averages of u 2 n as [M n (t 2 ) — M n (ti)}/(t 2 — ii), where 

dM n 



dt 



Mr 



(B.4) 
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Energy conservation implies that Y.n=o M n should grow linearly with time; 
hence we must use conservative values of u n on the right-hand side of B.4. 
Since we only compute conservative solutions u n at the beginning and the end 
of each time step, we are restricted to using a (second-order) trapezoidal rule 
to perform the integration of (B.4). 
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